rm(list=ls())

#---------------------------------------------------------------#
# PROJECT: HUMAN RIGHTS ABUSES                                  #
# SCRIPT: REPLICATION MANUSCRIPT - FIG 2,3,4,7,8,9              #
# AUTHOR: FLORES MACIAS AND ZARKIN                              #
# DATE: JANUARY  2022                                           #
#---------------------------------------------------------------#

#### ------ STEP 1: LOAD PACKAGES #### 

library(foreign)
library(tidyverse)
library(readstata13)
library(ggplot2)
library(PanelMatch)

#### ------ STEP 2: DEFINE THEME FOR FIGURES #### 


theme_bw1 <- function(base_size = 13, base_family = "") {
  theme_grey(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = 8, colour = "black",  hjust = .5 , vjust=1),
      axis.text.y =       element_text(size = 11 , colour = "black", hjust = 0 , vjust=.5 ), # changes position of X axis text
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90),
      axis.title.x =      element_text(size =12),
      legend.position = "bottom",
      panel.background = element_rect(fill = "white"),
      strip.background = element_rect(colour="black", fill="white")
    )
}

#### ------ STEP 3: LOAD DATA SETS #### 

hr<- read.dta13("0_DATA/ForMatching.dta")
hra<- read.dta13("0_DATA/ForMatching_Always.dta")
hrd<- read.dta13("0_DATA/ForMatching_NatDisOnly.dta")
hrd <- subset(hrd, year >= 2010)

#### ------ STEP 4: I, T IN INTEGER FORM #### 

hr$year <- as.integer(hr$year)
hr$inegi <- as.integer(hr$inegi)

hra$year <- as.integer(hra$year)
hra$inegi <- as.integer(hra$inegi)

hrd$year <- as.integer(hrd$year)
hrd$inegi <- as.integer(hrd$inegi)


#### ------ STEP 5: SUPPLEMENTAL APPENDIX FIGURE 3A ####

trfignd<- DisplayTreatment(unit.id = "inegi",
                           time.id = "year", legend.position = "bottom",
                           xlab = "Year", ylab = "Municipality",
                           treatment = "disaster_wo_operativo", data = hrd,
                           color.of.treated = "#636363",
                           color.of.untreated = "#f0f0f0",
                           legend.labels = c("Not Treated", "Treated"),
                           hide.y.tick.label = TRUE,
                           dense.plot = TRUE,
                           x.angle=90)
trfignd<-trfignd + labs(fill="", title="")
plot(trfignd)

ggsave("2_FIGURES/SA_FIG3A.png", plot=trfignd, width=10,height=8, units = "in", dpi = 300)


#### ------ STEP 6: SUPPLEMENTAL APPENDIX FIGURE 4A ####
##------------ DATA ####
msetsopps <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                        treatment = "operativo", refinement.method = "ps.weight", 
                        data = hr, match.missing = T, 
                        covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                        qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                        exact.match.variables = c("gradorezagosoc2010num", "rural"),
                        lead = 0:3, forbid.treatment.reversal = FALSE)

plot(msetsopps$att)
summary(msetsopps$att)

msopps<-data.frame(unlist(msetsopps))
msopps<-tibble::rownames_to_column(msopps, "unit") 
msopps <- data.frame(do.call('rbind', strsplit(as.character(msopps$unit),'.',fixed=TRUE)))
msopps <- transform(msopps, year = substr(X3, 1, 4), count = substr(X3, 5, 8))
msopps$inegiyear <- with(msopps, interaction(X2,  year))      
msopps <-aggregate(msopps, by = list(unique.values = msopps$inegiyear), 
                   FUN = length)
msopps$count2<-1

msopps2 <- msopps %>% group_by(inegiyear) %>%
  summarize(sum_count = sum(count2))

msopps2<-msopps2 %>% 
  rename(
    Matches=inegiyear,
    Observations=sum_count
  )



##------------ FIGURE ####

mset_figop<-ggplot(data=msopps2, aes(x=Matches, y=Observations)) +
  geom_bar(stat="identity", width=1) +
  ylab("Frequency of size") + 
  xlab("Matched set size") +
  theme_bw1() +
  xlim(0, 500)

plot(mset_figop)

ggsave("2_FIGURES/SA_FIG4A.png", plot=mset_figop, width=7,height=4, units = "in", dpi = 300)

rm(msetsopps,msopps)

#### ------ STEP 7: SUPPLEMENTAL APPENDIX FIGURE 5A ####
##------------ DATA ####
natdis_match <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                           treatment = "disaster_wo_operativo", refinement.method = "ps.weight", 
                           data = hrd, match.missing = T, 
                           covs.formula = ~ I(lag(lnpob, 1:2)) + I(lag(subsetSECFORCESpc, 1:2)) , 
                           qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                           exact.match.variables = c("gradorezagosoc2010num", "rural"),
                           lead = 0:2, forbid.treatment.reversal = FALSE)

plot(natdis_match$att)
summary(natdis_match$att)

msnatdis<-data.frame(unlist(natdis_match))
msnatdis<-tibble::rownames_to_column(msnatdis, "unit") 
msnatdis <- data.frame(do.call('rbind', strsplit(as.character(msnatdis$unit),'.',fixed=TRUE)))
msnatdis <- transform(msnatdis, year = substr(X3, 1, 4), count = substr(X3, 5, 8))
msnatdis$inegiyear <- with(msnatdis, interaction(X2,  year))      
msnatdis <-aggregate(msnatdis, by = list(unique.values = msnatdis$inegiyear), 
                     FUN = length)
msnatdis$count2<-1

msnatdis2 <- msnatdis %>% group_by(inegiyear) %>%
  summarize(sum_count = sum(count2))

msnatdis2<-msnatdis2 %>% 
  rename(
    Matches=inegiyear,
    Observations=sum_count
  )

##------------ FIGURE ####

mset_fignd<-ggplot(data=msnatdis2, aes(x=Matches, y=Observations)) +
  geom_bar(stat="identity", width=1) +
  ylab("Frequency of size") + 
  xlab("Matched set size") +
  theme_bw1() +
  xlim(0, 300)

plot(mset_fignd)

ggsave("2_FIGURES/SA_FIG5A.png", plot=mset_fignd, width=7,height=4, units = "in", dpi = 300)

rm(msnatdis, natdis_match)

#### ------ STEP 8: SUPPLEMENTAL APPENDIX FIGURE 6A #### 
##------------ PS WEIGHTS ####


# MATCHING
otherps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                            treatment = "operativo", refinement.method = "ps.weight", 
                            data = hr, match.missing = T, 
                            covs.formula = ~ I(lag(quejasCNDHpc, 1:3))  + I(lag(pob15a29pc, 1:3)) + I(lag(changeparty, 1:3)) + I(lag(isal, 1:3)) + I(lag(ieduc, 1:3)) + I(lag(iing, 1:3)), 
                            exact.match.variables = c("gradorezagosoc2010num", "rural"),
                            qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                            lead = 0:3, forbid.treatment.reversal = FALSE)


# CREATE DATA SET
cov_otherps<-data.frame(get_covariate_balance(otherps_match$att, 
                                              hr, covariates = c("quejasCNDHpc", "pob15a29pc", 
                                                                 "changeparty", "isal", "ieduc", "iing"), 
                                              plot = F, ylim = c(-2,2)))

##------------ CBPS WEIGHTS ####


othercbps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                              treatment = "operativo", refinement.method = "CBPS.weight", 
                              data = hr, match.missing = T, 
                              covs.formula = ~ I(lag(quejasCNDHpc, 1:3))  + I(lag(pob15a29pc, 1:3)) + I(lag(changeparty, 1:3)) + I(lag(isal, 1:3)) + I(lag(ieduc, 1:3)) + I(lag(iing, 1:3)), 
                              qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                              exact.match.variables = c("gradorezagosoc2010num", "rural"),
                              lead = 0:3, forbid.treatment.reversal = FALSE)

# CREATE DATA SET
cov_othercbps<-data.frame(get_covariate_balance(othercbps_match$att, 
                                                hr, covariates = c("quejasCNDHpc", "pob15a29pc", 
                                                                   "changeparty", "isal", "ieduc", "iing"), 
                                                plot = F, ylim = c(-2,2)))

##------------ FIGURE ####


# GENERATE NEW VARIABLES AND RENAME

cov_othercbps$Time<-c("1","2","3")
cov_otherps$Time<-c("1","2","3")


cov_othercbps<-cov_othercbps %>% 
  pivot_longer(!Time, names_to = "Covariate", values_to = "Estimate")
cov_otherps<-cov_otherps %>% 
  pivot_longer(!Time, names_to = "Covariate", values_to = "Estimate")

cov_othercbps$Method<-c("CBPS Weights")
cov_otherps$Method<-c("PS Weights")


cov_othercbps$Covariate<-recode(cov_othercbps$Covariate, "quejasCNDHpc" = "All complaints rate", "pob15a29pc" = "% 15-to-29 year olds", "changeparty" = "Change party mayor", 
                                "isal" = "Health Index", "ieduc" = "Education Index", "iing" = "Income Index")
cov_otherps$Covariate<-recode(cov_otherps$Covariate, "quejasCNDHpc" = "All complaints rate", "pob15a29pc" = "% 15-to-29 year olds", "changeparty" = "Change party mayor", 
                              "isal" = "Health Index", "ieduc" = "Education Index", "iing" = "Income Index")

# APPEND DATA SETS

cov_other<-rbind(cov_othercbps,cov_otherps)

cov_other$Time<-as.numeric(cov_other$Time)

# CREATE FIGURE

cov_otherfig=ggplot(cov_other, aes(x=Time, y=Estimate, color=Covariate)) +
  geom_line(size=1) +
  facet_wrap(vars(Method)) +
  theme_bw1() +
  geom_hline(yintercept=0, color="black",  linetype=2) +
  scale_x_continuous(breaks=c(1,2,3), labels=c("t-3", "t-2", "t-1")) +
  scale_y_continuous(name="Mean difference (st. dev)", limits=c(-2,2)) +
  scale_color_grey()

plot(cov_otherfig)

ggsave("2_FIGURES/SA_FIG6A.png", plot=cov_otherfig, width=10,height=5, units = "in", dpi = 300)

rm(cov_otherps, cov_othercbps, othercbps_match, otherps_match)

### ------ STEP 9: SUPPLEMENTAL APPENDIX FIGURE 7A #### 
##------------ NO MATCHING METHOD ####

# MATCHING
ndmatch_none <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                           treatment = "disaster_wo_operativo", refinement.method = "none", 
                           data = hrd, match.missing = T, 
                           covs.formula = ~ I(lag(lnpob, 1:2)) + I(lag(subsetSECFORCESpc, 1:2)) , 
                           qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                           exact.match.variables = c("gradorezagosoc2010num", "rural"),
                           lead = 0:2, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE

get_covariate_balance(ndmatch_none$att, hrd, covariates = c("lnpob", "subsetSECFORCESpc"), plot = F, ylim = c(-2,2))

# CREATE DATA SET
ndcov_none<-data.frame(get_covariate_balance(ndmatch_none$att, hrd, covariates = 
                                               c("lnpob", "subsetSECFORCESpc"), 
                                             plot = F, ylim = c(-2,2)))

##------------ PS MATCH 5 ####

# MATCHING
ndmatch_ps5 <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                          treatment = "disaster_wo_operativo", refinement.method = "ps.match", 
                          data = hrd, match.missing = T, 
                          covs.formula = ~ I(lag(lnpob, 1:2)) + I(lag(subsetSECFORCESpc, 1:2)) , 
                          size.match = 5, qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                          exact.match.variables = c("gradorezagosoc2010num", "rural"),
                          lead = 0:2, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE

get_covariate_balance(ndmatch_ps5$att, hrd, covariates = c("lnpob", "subsetSECFORCESpc"), plot = F, ylim = c(-2,2))

# CREATE DATA SET
ndcov_ps5<-data.frame(get_covariate_balance(ndmatch_ps5$att, hrd, covariates = 
                                              c("lnpob", "subsetSECFORCESpc"), 
                                            plot = F, ylim = c(-2,2)))


##------------ PS MATCH 10 ####

# MATCHING
ndmatch_ps10 <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                           treatment = "disaster_wo_operativo", refinement.method = "ps.match", 
                           data = hrd, match.missing = T, 
                           covs.formula = ~ I(lag(lnpob, 1:2)) + I(lag(subsetSECFORCESpc, 1:2)) , 
                           size.match = 10, qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                           exact.match.variables = c("gradorezagosoc2010num", "rural"),
                           lead = 0:2, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE

get_covariate_balance(ndmatch_ps10$att, hrd, covariates = c("lnpob", "subsetSECFORCESpc"), plot = F, ylim = c(-2,2))

# CREATE DATA SET
ndcov_ps10<-data.frame(get_covariate_balance(ndmatch_ps10$att, hrd, covariates = 
                                               c("lnpob", "subsetSECFORCESpc"), 
                                             plot = F, ylim = c(-2,2)))

##------------ CBPS MATCH 5 ####

# MATCHING
ndmatch_cbps5 <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                            treatment = "disaster_wo_operativo", refinement.method = "CBPS.match", 
                            data = hrd, match.missing = T, 
                            covs.formula = ~ I(lag(lnpob, 1:2)) + I(lag(subsetSECFORCESpc, 1:2)) , 
                            size.match = 5, qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                            exact.match.variables = c("gradorezagosoc2010num", "rural"),
                            lead = 0:2, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE

get_covariate_balance(ndmatch_cbps5$att, hrd, covariates = c("lnpob", "subsetSECFORCESpc"), plot = F, ylim = c(-2,2))

# CREATE DATA SET
ndcov_cbps5<-data.frame(get_covariate_balance(ndmatch_cbps5$att, hrd, covariates = 
                                                c("lnpob", "subsetSECFORCESpc"), 
                                              plot = F, ylim = c(-2,2)))

##------------ CBPS MATCH 10 ####

# MATCHING
ndmatch_cbps10 <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                             treatment = "disaster_wo_operativo", refinement.method = "CBPS.match", 
                             data = hrd, match.missing = T, 
                             covs.formula = ~ I(lag(lnpob, 1:2)) + I(lag(subsetSECFORCESpc, 1:2)) , 
                             size.match = 10, qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                             exact.match.variables = c("gradorezagosoc2010num", "rural"),
                             lead = 0:2, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE

get_covariate_balance(ndmatch_cbps10$att, hrd, covariates = c("lnpob", "subsetSECFORCESpc"), plot = F, ylim = c(-2,2))

# CREATE DATA SET
ndcov_cbps10<-data.frame(get_covariate_balance(ndmatch_cbps10$att, hrd, covariates = 
                                                 c("lnpob", "subsetSECFORCESpc"), 
                                               plot = F, ylim = c(-2,2)))

##------------ MAHALANOBIS ####

# MATCHING
ndmatch_mahal <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                            treatment = "disaster_wo_operativo", refinement.method = "mahalanobis", 
                            data = hrd, match.missing = T, 
                            covs.formula = ~ I(lag(lnpob, 1:2)) + I(lag(subsetSECFORCESpc, 1:2)) , 
                            qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                            exact.match.variables = c("gradorezagosoc2010num", "rural"),
                            lead = 0:2, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE

get_covariate_balance(ndmatch_mahal$att, hrd, covariates = c("lnpob", "subsetSECFORCESpc"), plot = F, ylim = c(-2,2))

# CREATE DATA SET
ndcov_mahal<-data.frame(get_covariate_balance(ndmatch_mahal$att, hrd, covariates = 
                                                c("lnpob", "subsetSECFORCESpc"), 
                                              plot = F, ylim = c(-2,2)))

##------------ PS WEIGHTS ####

# MATCHING
ndmatch_psw <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                          treatment = "disaster_wo_operativo", refinement.method = "ps.weight", 
                          data = hrd, match.missing = T, 
                          covs.formula = ~ I(lag(lnpob, 1:2)) + I(lag(subsetSECFORCESpc, 1:2)) , 
                          qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                          exact.match.variables = c("gradorezagosoc2010num", "rural"),
                          lead = 0:2, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE

get_covariate_balance(ndmatch_psw$att, hrd, covariates = c("lnpob", "subsetSECFORCESpc"), plot = F, ylim = c(-2,2))

# CREATE DATA SET
ndcov_psw<-data.frame(get_covariate_balance(ndmatch_psw$att, hrd, covariates = 
                                              c("lnpob", "subsetSECFORCESpc"), 
                                            plot = F, ylim = c(-2,2)))

##------------ CBPS WEIGHTS ####

# MATCHING
ndmatch_cbpsw <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                            treatment = "disaster_wo_operativo", refinement.method = "CBPS.weight", 
                            data = hrd, match.missing = T, 
                            covs.formula = ~ I(lag(lnpob, 1:2)) + I(lag(subsetSECFORCESpc, 1:2)) , 
                            qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                            exact.match.variables = c("gradorezagosoc2010num", "rural"),
                            lead = 0:2, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE

get_covariate_balance(ndmatch_cbpsw$att, hrd, covariates = c("lnpob", "subsetSECFORCESpc"), plot = F, ylim = c(-2,2))

# CREATE DATA SET
ndcov_cbpsw<-data.frame(get_covariate_balance(ndmatch_cbpsw$att, hrd, covariates = 
                                                c("lnpob", "subsetSECFORCESpc"), 
                                              plot = F, ylim = c(-2,2)))

##------------ FIGURE ####


# GENERATE NEW VARIABLES AND RENAME

ndcov_cbps10$Time<-c("1","2")
ndcov_cbps5$Time<-c("1","2")
ndcov_cbpsw$Time<-c("1","2")
ndcov_mahal$Time<-c("1","2")
ndcov_none$Time<-c("1","2")
ndcov_ps10$Time<-c("1","2")
ndcov_ps5$Time<-c("1","2")
ndcov_psw$Time<-c("1","2")

ndcov_cbps10<-ndcov_cbps10 %>% 
  pivot_longer(!Time, names_to = "Covariate", values_to = "Estimate")
ndcov_cbps5<-ndcov_cbps5 %>% 
  pivot_longer(!Time, names_to = "Covariate", values_to = "Estimate")
ndcov_cbpsw<-ndcov_cbpsw %>% 
  pivot_longer(!Time, names_to = "Covariate", values_to = "Estimate")
ndcov_mahal<-ndcov_mahal %>% 
  pivot_longer(!Time, names_to = "Covariate", values_to = "Estimate")
ndcov_none<-ndcov_none %>% 
  pivot_longer(!Time, names_to = "Covariate", values_to = "Estimate")
ndcov_ps10<-ndcov_ps10 %>% 
  pivot_longer(!Time, names_to = "Covariate", values_to = "Estimate")
ndcov_ps5<-ndcov_ps5 %>% 
  pivot_longer(!Time, names_to = "Covariate", values_to = "Estimate")
ndcov_psw<-ndcov_psw %>% 
  pivot_longer(!Time, names_to = "Covariate", values_to = "Estimate")



ndcov_cbps10$Method<-c("CBPS match 10")
ndcov_cbps5$Method<-c("CBPS match 5")
ndcov_cbpsw$Method<-c("CBPS weights")
ndcov_mahal$Method<-c("Mahalanobis")
ndcov_none$Method<-c("None")
ndcov_ps10$Method<-c("PS match 10")
ndcov_ps5$Method<-c("PS match 5")
ndcov_psw$Method<-c("PS weights")


ndcov_cbps10$Covariate<-recode(ndcov_cbps10$Covariate, "subsetSECFORCESpc" = "Complaints rate", "lnpob" = "Log population")
ndcov_cbps5$Covariate<-recode(ndcov_cbps5$Covariate,"subsetSECFORCESpc" = "Complaints rate", "lnpob" = "Log population")
ndcov_cbpsw$Covariate<-recode(ndcov_cbpsw$Covariate,"subsetSECFORCESpc" = "Complaints rate", "lnpob" = "Log population")
ndcov_mahal$Covariate<-recode(ndcov_mahal$Covariate,"subsetSECFORCESpc" = "Complaints rate", "lnpob" = "Log population")
ndcov_none$Covariate<-recode(ndcov_none$Covariate,"subsetSECFORCESpc" = "Complaints rate", "lnpob" = "Log population")
ndcov_ps10$Covariate<-recode(ndcov_ps10$Covariate,"subsetSECFORCESpc" = "Complaints rate", "lnpob" = "Log population")
ndcov_ps5$Covariate<-recode(ndcov_ps5$Covariate,"subsetSECFORCESpc" = "Complaints rate", "lnpob" = "Log population")
ndcov_psw$Covariate<-recode(ndcov_psw$Covariate,"subsetSECFORCESpc" = "Complaints rate", "lnpob" = "Log population")


# APPEND DATA SETS

cov_nd<-rbind(ndcov_cbps10,ndcov_cbps5,ndcov_cbpsw,ndcov_mahal,
              ndcov_none, ndcov_ps10, ndcov_ps5, ndcov_psw)

cov_nd$Method <- factor(cov_nd$Method, levels = c("None", "PS weights", "PS match 5", "PS match 10", 
                                                  
                                                  "CBPS weights", "CBPS match 5", "CBPS match 10", "Mahalanobis"))
cov_nd$Time<-as.numeric(cov_nd$Time)

# CREATE FIGURE

covndfig=ggplot(cov_nd, aes(x=Time, y=Estimate, color=Covariate)) +
  geom_line(size=1) +
  facet_wrap(vars(Method), nrow=2) +
  theme_bw1() +
  geom_hline(yintercept=0, color="black",  linetype=2) +
  scale_x_continuous(breaks=c(1,2), labels=c("t-2", "t-1")) +
  scale_y_continuous(name="Mean difference (st. dev)", limits=c(-2,2)) +
  scale_color_grey()

plot(covndfig)

ggsave("2_FIGURES/SA_FIG7A.png", plot=covndfig, width=12,height=6, units = "in", dpi = 300)

rm(ndcov_cbps10,ndcov_cbps5,ndcov_cbpsw,ndcov_mahal,ndcov_none,ndcov_ps10,ndcov_ps5,
   ndcov_psw,ndmatch_cbps10,ndmatch_cbps5,ndmatch_cbpsw,ndmatch_mahal,ndmatch_none,
   ndmatch_ps10,ndmatch_ps5,ndmatch_psw)

#### ------ STEP 10: SUPPLEMENTAL APPENDIX FIGURE 8A ####
##------------ PS WEIGHTS ####

# 1. ONE LAG, ZERO LEAD

oneps_match <- PanelMatch(lag = 1, time.id = "year", unit.id = "inegi", 
                          treatment = "operativo", refinement.method = "ps.weight", 
                          data = hr, match.missing = T, 
                          covs.formula = ~ I(lag(thom, 1:1))  + I(lag(subsetSECFORCESpc, 1:1)) + I(lag(lnpob, 1:1)) + I(lag(turfwar_threeaf, 1:1)) + I(lag(samegovpres, 1:1)), 
                          qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                          exact.match.variables = c("gradorezagosoc2010num", "rural"),
                          lead = 0:0, forbid.treatment.reversal = FALSE)

oneps_estim <- PanelEstimate(sets = oneps_match, data = hr)

summary(oneps_estim)
one_ps<-as.data.frame(summary(oneps_estim))

# 2. ONE LAG, ONE LEADS

twops_match <- PanelMatch(lag = 1, time.id = "year", unit.id = "inegi", 
                          treatment = "operativo", refinement.method = "ps.weight", 
                          data = hr, match.missing = T, 
                          covs.formula = ~ I(lag(thom, 1:1))  + I(lag(subsetSECFORCESpc, 1:1)) + I(lag(lnpob, 1:1)) + I(lag(turfwar_threeaf, 1:1)) + I(lag(samegovpres, 1:1)), 
                          qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                          exact.match.variables = c("gradorezagosoc2010num", "rural"),
                          lead = 0:1, forbid.treatment.reversal = FALSE)

twops_estim <- PanelEstimate(sets = twops_match, data = hr)

summary(twops_estim)
two_ps<-as.data.frame(summary(twops_estim))

# 3. ONE LAG, TWO LEADS

threeps_match <- PanelMatch(lag = 1, time.id = "year", unit.id = "inegi", 
                            treatment = "operativo", refinement.method = "ps.weight", 
                            data = hr, match.missing = T, 
                            covs.formula = ~ I(lag(thom, 1:1))  + I(lag(subsetSECFORCESpc, 1:1)) + I(lag(lnpob, 1:1)) + I(lag(turfwar_threeaf, 1:1)) + I(lag(samegovpres, 1:1)), 
                            qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                            exact.match.variables = c("gradorezagosoc2010num", "rural"),
                            lead = 0:2, forbid.treatment.reversal = FALSE)

threeps_estim <- PanelEstimate(sets = threeps_match, data = hr)

summary(threeps_estim)
three_ps<-as.data.frame(summary(threeps_estim))

# 4. ONE LAG, THREE LEADS

fourps_match <- PanelMatch(lag = 1, time.id = "year", unit.id = "inegi", 
                           treatment = "operativo", refinement.method = "ps.weight", 
                           data = hr, match.missing = T, 
                           covs.formula = ~ I(lag(thom, 1:1))  + I(lag(subsetSECFORCESpc, 1:1)) + I(lag(lnpob, 1:1)) + I(lag(turfwar_threeaf, 1:1)) + I(lag(samegovpres, 1:1)), 
                           qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                           exact.match.variables = c("gradorezagosoc2010num", "rural"),
                           lead = 0:3, forbid.treatment.reversal = FALSE)

fourps_estim <- PanelEstimate(sets = fourps_match, data = hr)

summary(fourps_estim)
four_ps<-as.data.frame(summary(fourps_estim))


# 5. TWO LAG, ZERO LEAD

fiveps_match <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                           treatment = "operativo", refinement.method = "ps.weight", 
                           data = hr, match.missing = T, 
                           covs.formula = ~ I(lag(thom, 1:2))  + I(lag(subsetSECFORCESpc, 1:2)) + I(lag(lnpob, 1:2)) + I(lag(turfwar_threeaf, 1:2)) + I(lag(samegovpres, 1:2)), 
                           qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                           exact.match.variables = c("gradorezagosoc2010num", "rural"),
                           lead = 0:0, forbid.treatment.reversal = FALSE)

fiveps_estim <- PanelEstimate(sets = fiveps_match, data = hr)

summary(fiveps_estim)
five_ps<-as.data.frame(summary(fiveps_estim))


# 6. TWO LAG, ONE LEAD

sixps_match <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                          treatment = "operativo", refinement.method = "ps.weight", 
                          data = hr, match.missing = T, 
                          covs.formula = ~ I(lag(thom, 1:2))  + I(lag(subsetSECFORCESpc, 1:2)) + I(lag(lnpob, 1:2)) + I(lag(turfwar_threeaf, 1:2)) + I(lag(samegovpres, 1:2)), 
                          qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                          exact.match.variables = c("gradorezagosoc2010num", "rural"),
                          lead = 0:1, forbid.treatment.reversal = FALSE)

sixps_estim <- PanelEstimate(sets = sixps_match, data = hr)

summary(sixps_estim)
six_ps<-as.data.frame(summary(sixps_estim))



# 7. TWO LAG, TWO LEAD

sevenps_match <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                            treatment = "operativo", refinement.method = "ps.weight", 
                            data = hr, match.missing = T, 
                            covs.formula = ~ I(lag(thom, 1:2))  + I(lag(subsetSECFORCESpc, 1:2)) + I(lag(lnpob, 1:2)) + I(lag(turfwar_threeaf, 1:2)) + I(lag(samegovpres, 1:2)), 
                            qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                            exact.match.variables = c("gradorezagosoc2010num", "rural"),
                            lead = 0:2, forbid.treatment.reversal = FALSE)

sevenps_estim <- PanelEstimate(sets = sevenps_match, data = hr)

summary(sevenps_estim)
seven_ps<-as.data.frame(summary(sevenps_estim))


# 8. TWO LAG, THREE LEAD

eightps_match <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                            treatment = "operativo", refinement.method = "ps.weight", 
                            data = hr, match.missing = T, 
                            covs.formula = ~ I(lag(thom, 1:2))  + I(lag(subsetSECFORCESpc, 1:2)) + I(lag(lnpob, 1:2)) + I(lag(turfwar_threeaf, 1:2)) + I(lag(samegovpres, 1:2)), 
                            qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                            exact.match.variables = c("gradorezagosoc2010num", "rural"),
                            lead = 0:3, forbid.treatment.reversal = FALSE)

eightps_estim <- PanelEstimate(sets = eightps_match, data = hr)

summary(eightps_estim)
eight_ps<-as.data.frame(summary(eightps_estim))


# 9. THREE LAG, ZERO LEAD

nineps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                           treatment = "operativo", refinement.method = "ps.weight", 
                           data = hr, match.missing = T, 
                           covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                           qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                           exact.match.variables = c("gradorezagosoc2010num", "rural"),
                           lead = 0:0, forbid.treatment.reversal = FALSE)

nineps_estim <- PanelEstimate(sets = nineps_match, 
                              data = hr)
summary(nineps_estim)
nine_ps<-as.data.frame(summary(nineps_estim))


# 10. THREE LAG, ONE LEAD

tenps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                          treatment = "operativo", refinement.method = "ps.weight", 
                          data = hr, match.missing = T, 
                          covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                          qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                          exact.match.variables = c("gradorezagosoc2010num", "rural"),
                          lead = 0:1, forbid.treatment.reversal = FALSE)

tenps_estim <- PanelEstimate(sets = tenps_match, data = hr)

summary(tenps_estim)
ten_ps<-as.data.frame(summary(tenps_estim))

# 11. THREE LAG, TWO LEAD

elevenps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                             treatment = "operativo", refinement.method = "ps.weight", 
                             data = hr, match.missing = T, 
                             covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                             qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                             exact.match.variables = c("gradorezagosoc2010num", "rural"),
                             lead = 0:2, forbid.treatment.reversal = FALSE)

elevenps_estim <- PanelEstimate(sets = elevenps_match, data = hr)

summary(elevenps_estim)
eleven_ps<-as.data.frame(summary(elevenps_estim))


# 12. THREE LAG, THREE LEAD

twelveps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                             treatment = "operativo", refinement.method = "ps.weight", 
                             data = hr, match.missing = T, 
                             covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                             qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                             exact.match.variables = c("gradorezagosoc2010num", "rural"),
                             lead = 0:3, forbid.treatment.reversal = FALSE)

twelveps_estim <- PanelEstimate(sets = twelveps_match, data = hr)

summary(twelveps_estim)
twelve_ps<-as.data.frame(summary(twelveps_estim))

rm(eightps_estim,eightps_match,fiveps_estim,fiveps_match,fourps_estim,fourps_match,
   nineps_estim,nineps_match,oneps_estim,oneps_match,sevenps_estim,sevenps_match,sixps_estim,sixps_match,
   tenps_estim,tenps_match,threeps_estim,threeps_match,twelveps_estim,twelveps_match,elevenps_estim,elevenps_match,
   twops_estim,twops_match)

##------------ CBPS WEIGHTS ####

# 1. ONE LAG, ZERO LEAD

onecpbs_match <- PanelMatch(lag = 1, time.id = "year", unit.id = "inegi", 
                            treatment = "operativo", refinement.method = "CBPS.weight", 
                            data = hr, match.missing = T, 
                            covs.formula = ~ I(lag(thom, 1:1))  + I(lag(subsetSECFORCESpc, 1:1)) + I(lag(lnpob, 1:1)) + I(lag(turfwar_threeaf, 1:1)) + I(lag(samegovpres, 1:1)), 
                            size.match = 5, qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                            exact.match.variables = c("gradorezagosoc2010num", "rural"),
                            lead = 0:0, forbid.treatment.reversal = FALSE)

onecbps_estim <- PanelEstimate(sets = onecpbs_match, data = hr)

summary(onecbps_estim)
one_cbps<-as.data.frame(summary(onecbps_estim))

# 2. ONE LAG, ONE LEADS

twocbps_match <- PanelMatch(lag = 1, time.id = "year", unit.id = "inegi", 
                            treatment = "operativo", refinement.method = "CBPS.weight", 
                            data = hr, match.missing = T, 
                            covs.formula = ~ I(lag(thom, 1:1))  + I(lag(subsetSECFORCESpc, 1:1)) + I(lag(lnpob, 1:1)) + I(lag(turfwar_threeaf, 1:1)) + I(lag(samegovpres, 1:1)), 
                            size.match = 5, qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                            exact.match.variables = c("gradorezagosoc2010num", "rural"),
                            lead = 0:1, forbid.treatment.reversal = FALSE)

twocbps_estim <- PanelEstimate(sets = twocbps_match, data = hr)

summary(twocbps_estim)
two_cbps<-as.data.frame(summary(twocbps_estim)) 

# 3. ONE LAG, TWO LEADS

threecbps_match <- PanelMatch(lag = 1, time.id = "year", unit.id = "inegi", 
                              treatment = "operativo", refinement.method = "CBPS.weight", 
                              data = hr, match.missing = T, 
                              covs.formula = ~ I(lag(thom, 1:1))  + I(lag(subsetSECFORCESpc, 1:1)) + I(lag(lnpob, 1:1)) + I(lag(turfwar_threeaf, 1:1)) + I(lag(samegovpres, 1:1)), 
                              qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                              exact.match.variables = c("gradorezagosoc2010num", "rural"),
                              lead = 0:2, forbid.treatment.reversal = FALSE)

threecbps_estim <- PanelEstimate(sets = threecbps_match, data = hr)

summary(threecbps_estim)
three_cbps<-as.data.frame(summary(threecbps_estim))

# 4. ONE LAG, THREE LEADS

fourcbps_match <- PanelMatch(lag = 1, time.id = "year", unit.id = "inegi", 
                             treatment = "operativo", refinement.method = "CBPS.weight", 
                             data = hr, match.missing = T, 
                             covs.formula = ~ I(lag(thom, 1:1))  + I(lag(subsetSECFORCESpc, 1:1)) + I(lag(lnpob, 1:1)) + I(lag(turfwar_threeaf, 1:1)) + I(lag(samegovpres, 1:1)), 
                             qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                             exact.match.variables = c("gradorezagosoc2010num", "rural"),
                             lead = 0:3, forbid.treatment.reversal = FALSE)

fourcbps_estim <- PanelEstimate(sets = fourcbps_match, data = hr)

summary(fourcbps_estim)
four_cbps<-as.data.frame(summary(fourcbps_estim))

# 5. TWO LAG, ZERO LEAD

fivecbps_match <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                             treatment = "operativo", refinement.method = "CBPS.weight", 
                             data = hr, match.missing = T, 
                             covs.formula = ~ I(lag(thom, 1:2))  + I(lag(subsetSECFORCESpc, 1:2)) + I(lag(lnpob, 1:2)) + I(lag(turfwar_threeaf, 1:2)) + I(lag(samegovpres, 1:2)), 
                             qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                             exact.match.variables = c("gradorezagosoc2010num", "rural"),
                             lead = 0:0, forbid.treatment.reversal = FALSE)

fivecbps_estim <- PanelEstimate(sets = fivecbps_match, data = hr)

summary(fivecbps_estim)
five_cbps<-as.data.frame(summary(fivecbps_estim))

# 6. TWO LAG, ONE LEAD

sixcbps_match <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                            treatment = "operativo", refinement.method = "CBPS.weight", 
                            data = hr, match.missing = T, 
                            covs.formula = ~ I(lag(thom, 1:2))  + I(lag(subsetSECFORCESpc, 1:2)) + I(lag(lnpob, 1:2)) + I(lag(turfwar_threeaf, 1:2)) + I(lag(samegovpres, 1:2)), 
                            qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                            exact.match.variables = c("gradorezagosoc2010num", "rural"),
                            lead = 0:1, forbid.treatment.reversal = FALSE)

sixcbps_estim <- PanelEstimate(sets = sixcbps_match, data = hr)

summary(sixcbps_estim)
six_cbps<-as.data.frame(summary(sixcbps_estim))


# 7. TWO LAG, TWO LEAD

sevencbps_match <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                              treatment = "operativo", refinement.method = "CBPS.weight", 
                              data = hr, match.missing = T, 
                              covs.formula = ~ I(lag(thom, 1:2))  + I(lag(subsetSECFORCESpc, 1:2)) + I(lag(lnpob, 1:2)) + I(lag(turfwar_threeaf, 1:2)) + I(lag(samegovpres, 1:2)), 
                              qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                              exact.match.variables = c("gradorezagosoc2010num", "rural"),
                              lead = 0:2, forbid.treatment.reversal = FALSE)

sevencbps_estim <- PanelEstimate(sets = sevencbps_match, data = hr)

summary(sevencbps_estim)
seven_cbps<-as.data.frame(summary(sevencbps_estim))

# 8. TWO LAG, THREE LEAD

eightcbps_match <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                              treatment = "operativo", refinement.method = "CBPS.weight", 
                              data = hr, match.missing = T, 
                              covs.formula = ~ I(lag(thom, 1:2))  + I(lag(subsetSECFORCESpc, 1:2)) + I(lag(lnpob, 1:2)) + I(lag(turfwar_threeaf, 1:2)) + I(lag(samegovpres, 1:2)), 
                              qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                              exact.match.variables = c("gradorezagosoc2010num", "rural"),
                              lead = 0:3, forbid.treatment.reversal = FALSE)

eightcbps_estim <- PanelEstimate(sets = eightcbps_match, data = hr)

summary(eightcbps_estim)
eight_cbps<-as.data.frame(summary(eightcbps_estim))

# 9. THREE LAG, ZERO LEAD

ninecbps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                             treatment = "operativo", refinement.method = "CBPS.weight", 
                             data = hr, match.missing = T, 
                             covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                             qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                             exact.match.variables = c("gradorezagosoc2010num", "rural"),
                             lead = 0:0, forbid.treatment.reversal = FALSE)

ninecbps_estim <- PanelEstimate(sets = ninecbps_match, data = hr)

summary(ninecbps_estim)
nine_cbps<-as.data.frame(summary(ninecbps_estim))

# 10. THREE LAG, ONE LEAD

tencbps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                            treatment = "operativo", refinement.method = "CBPS.weight", 
                            data = hr, match.missing = T, 
                            covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                            qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                            exact.match.variables = c("gradorezagosoc2010num", "rural"),
                            lead = 0:1, forbid.treatment.reversal = FALSE)

tencbps_estim <- PanelEstimate(sets = tencbps_match, data = hr)

summary(tencbps_estim)
ten_cbps<-as.data.frame(summary(tencbps_estim))

# 11. THREE LAG, TWO LEAD

elevencbps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                               treatment = "operativo", refinement.method = "CBPS.weight", 
                               data = hr, match.missing = T, 
                               covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                               qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                               exact.match.variables = c("gradorezagosoc2010num", "rural"),
                               lead = 0:2, forbid.treatment.reversal = FALSE)

elevencbps_estim <- PanelEstimate(sets = elevencbps_match, data = hr)

summary(elevencbps_estim)
eleven_cbps<-as.data.frame(summary(elevencbps_estim))


# 12. THREE LAG, THREE LEAD

twelvecbps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                               treatment = "operativo", refinement.method = "CBPS.weight", 
                               data = hr, match.missing = T, 
                               covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                               qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                               exact.match.variables = c("gradorezagosoc2010num", "rural"),
                               lead = 0:3, forbid.treatment.reversal = FALSE)

twelvecbps_estim <- PanelEstimate(sets = twelvecbps_match, data = hr)

summary(twelvecbps_estim)
twelve_cbps<-as.data.frame(summary(twelvecbps_estim))

rm(onecpbs_match,onecbps_estim,twocbps_estim,twocbps_match,threecbps_estim,threecbps_match,fourcbps_estim,fourcbps_match,
   fivecbps_estim,fivecbps_match,sixcbps_estim,sixcbps_match,sevencbps_estim,sevencbps_match,
   eightcbps_estim,eightcbps_match,ninecbps_estim,ninecbps_match,tencbps_estim,tencbps_match,
   elevencbps_estim,elevencbps_match,twelvecbps_estim,twelvecbps_match)

##------------ FIGURE ####

# GENERATE NEW VARIABLES AND RENAME

one_ps$Method<-c("PS Weights")
one_cbps$Method<-c("CBPS Weights")
two_ps$Method<-c("PS Weights")
two_cbps$Method<-c("CBPS Weights")
three_ps$Method<-c("PS Weights")
three_cbps$Method<-c("CBPS Weights")
four_ps$Method<-c("PS Weights")
four_cbps$Method<-c("CBPS Weights")
five_ps$Method<-c("PS Weights")
five_cbps$Method<-c("CBPS Weights")
six_ps$Method<-c("PS Weights")
six_cbps$Method<-c("CBPS Weights")
seven_ps$Method<-c("PS Weights")
seven_cbps$Method<-c("CBPS Weights")
eight_ps$Method<-c("PS Weights")
eight_cbps$Method<-c("CBPS Weights")
nine_ps$Method<-c("PS Weights")
nine_cbps$Method<-c("CBPS Weights")
ten_ps$Method<-c("PS Weights")
ten_cbps$Method<-c("CBPS Weights")
eleven_ps$Method<-c("PS Weights")
eleven_cbps$Method<-c("CBPS Weights")
twelve_ps$Method<-c("PS Weights")
twelve_cbps$Method<-c("CBPS Weights")

one_ps<-one_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

one_cbps<-one_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

two_ps<-two_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

two_cbps<-two_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

three_ps<-three_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

three_cbps<-three_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

four_ps<-four_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

four_cbps<-four_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

five_ps<-five_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

five_cbps<-five_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

six_ps<-six_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

six_cbps<-six_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

seven_ps<-seven_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

seven_cbps<-seven_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

eight_ps<-eight_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

eight_cbps<-eight_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

nine_ps<-nine_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

nine_cbps<-nine_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

ten_ps<-ten_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

ten_cbps<-ten_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

eleven_ps<-eleven_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

eleven_cbps<-eleven_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

twelve_ps<-twelve_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

twelve_cbps<-twelve_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

one_ps$Time<-c("t+0")
one_cbps$Time<-c("t+0")
two_ps$Time<-c("t+0","t+1")
two_cbps$Time<-c("t+0","t+1")
three_ps$Time<-c("t+0","t+1","t+2")
three_cbps$Time<-c("t+0","t+1","t+2")
four_ps$Time<-c("t+0","t+1","t+2","t+3")
four_cbps$Time<-c("t+0","t+1","t+2","t+3")
five_ps$Time<-c("t+0")
five_cbps$Time<-c("t+0")
six_ps$Time<-c("t+0","t+1")
six_cbps$Time<-c("t+0","t+1")
seven_ps$Time<-c("t+0","t+1","t+2")
seven_cbps$Time<-c("t+0","t+1","t+2")
eight_ps$Time<-c("t+0","t+1","t+2","t+3")
eight_cbps$Time<-c("t+0","t+1","t+2","t+3")
nine_ps$Time<-c("t+0")
nine_cbps$Time<-c("t+0")
ten_ps$Time<-c("t+0","t+1")
ten_cbps$Time<-c("t+0","t+1")
eleven_ps$Time<-c("t+0","t+1","t+2")
eleven_cbps$Time<-c("t+0","t+1","t+2")
twelve_ps$Time<-c("t+0","t+1","t+2","t+3")
twelve_cbps$Time<-c("t+0","t+1","t+2","t+3")

one_ps$Model<-c("1 Lag 0 Lead")
one_cbps$Model<-c("1 Lag 0 Lead")
two_ps$Model<-c("1 Lag 1 Lead")
two_cbps$Model<-c("1 Lag 1 Lead")
three_ps$Model<-c("1 Lag 2 Lead")
three_cbps$Model<-c("1 Lag 2 Lead")
four_ps$Model<-c("1 Lag 3 Lead")
four_cbps$Model<-c("1 Lag 3 Lead")
five_ps$Model<-c("2 Lag 0 Lead")
five_cbps$Model<-c("2 Lag 0 Lead")
six_ps$Model<-c("2 Lag 1 Lead")
six_cbps$Model<-c("2 Lag 1 Lead")
seven_ps$Model<-c("2 Lag 2 Lead")
seven_cbps$Model<-c("2 Lag 2 Lead")
eight_ps$Model<-c("2 Lag 3 Lead")
eight_cbps$Model<-c("2 Lag 3 Lead")
nine_ps$Model<-c("3 Lag 0 Lead")
nine_cbps$Model<-c("3 Lag 0 Lead")
ten_ps$Model<-c("3 Lag 1 Lead")
ten_cbps$Model<-c("3 Lag 1 Lead")
eleven_ps$Model<-c("3 Lag 2 Lead")
eleven_cbps$Model<-c("3 Lag 2 Lead")
twelve_ps$Model<-c("3 Lag 3 Lead")
twelve_cbps$Model<-c("3 Lag 3 Lead")

# APPEND DATA SETS


lagleadops<-rbind(one_cbps,one_ps,two_cbps,two_ps,three_cbps,three_ps,
                  four_cbps,four_ps,five_cbps,five_ps,six_cbps,six_ps,
                  seven_cbps,seven_ps,eight_cbps,eight_ps,nine_cbps,nine_ps,
                  ten_cbps,ten_ps,eleven_cbps,eleven_ps,twelve_cbps,twelve_ps)

# CREATE FIGURE

theme_bw2 <- function(base_size = 13, base_family = "") {
  theme_grey(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = 8, colour = "black",  hjust = .5 , vjust=1),
      axis.text.y =       element_text(size = 11 , colour = "black", hjust = 0 , vjust=.5 ), # changes position of X axis text
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90),
      axis.title.x =      element_text(size =12),
      legend.position = "bottom",
      panel.background = element_rect(fill = "white"),
      strip.background = element_rect(colour="black", fill="white"),
      legend.title=element_blank(),
      legend.text=element_text(size=7),
      legend.key.size = unit(.7,"line")
    )
}

lagleadsops_fig=ggplot(lagleadops, aes(Time, Estimate, fill=Model)) + 
  geom_linerange(aes(ymin=CILow, ymax=CIHigh, colour=Model), size=1,  position = position_dodge(width = 0.5)) + 
  ylim(-2, 2) +
  theme_bw1() +
  geom_hline(yintercept=0, color="gray",  linetype=2) +
  facet_grid(.~Method)

plot(lagleadsops_fig)  

ggsave("2_FIGURES/SA_FIG8A.png", plot=lagleadsops_fig, width=8,height=4, units = "in", dpi = 300)

rm(one_cbps,one_ps,two_cbps,two_ps,three_cbps,three_ps,four_cbps,four_ps,
   five_cbps,five_ps,six_cbps,six_ps,seven_cbps,seven_ps,eight_cbps,eight_ps,
   nine_cbps,nine_ps,ten_cbps,ten_ps,eleven_cbps,eleven_ps,twelve_cbps,twelve_ps)

#### ------ STEP 11: SUPPLEMENTAL APPENDIX FIGURE 9A ####
##------------ PS WEIGHTS ####

# 1. ONE LAG, ZERO LEAD

onendps_match <- PanelMatch(lag = 1, time.id = "year", unit.id = "inegi", 
                            treatment = "disaster_wo_operativo", refinement.method = "ps.weight", 
                            data = hrd, match.missing = T, 
                            covs.formula = ~I(lag(subsetSECFORCESpc, 1:1)) + I(lag(lnpob, 1:1)) , 
                            qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                            exact.match.variables = c("gradorezagosoc2010num", "rural"),
                            lead = 0:0, forbid.treatment.reversal = FALSE)

onendps_estim <- PanelEstimate(sets = onendps_match, data = hrd)

summary(onendps_estim)
onend_ps<-as.data.frame(summary(onendps_estim))

# 2. ONE LAG, ONE LEADS

twondps_match <- PanelMatch(lag = 1, time.id = "year", unit.id = "inegi", 
                            treatment = "disaster_wo_operativo", refinement.method = "ps.weight", 
                            data = hrd, match.missing = T, 
                            covs.formula = ~ I(lag(subsetSECFORCESpc, 1:1)) + I(lag(lnpob, 1:1)) , 
                            qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                            exact.match.variables = c("gradorezagosoc2010num", "rural"),
                            lead = 0:1, forbid.treatment.reversal = FALSE)

twondps_estim <- PanelEstimate(sets = twondps_match, data = hrd)

summary(twondps_estim)
twond_ps<-as.data.frame(summary(twondps_estim))

# 3. ONE LAG, TWO LEADS

threendps_match <- PanelMatch(lag = 1, time.id = "year", unit.id = "inegi", 
                              treatment = "disaster_wo_operativo", refinement.method = "ps.weight", 
                              data = hrd, match.missing = T, 
                              covs.formula = ~ I(lag(subsetSECFORCESpc, 1:1)) + I(lag(lnpob, 1:1)) , 
                              qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                              exact.match.variables = c("gradorezagosoc2010num", "rural"),
                              lead = 0:2, forbid.treatment.reversal = FALSE)

threendps_estim <- PanelEstimate(sets = threendps_match, data = hrd)

summary(threendps_estim)
threend_ps<-as.data.frame(summary(threendps_estim))

# 4. TWO LAG, ZERO LEAD

fourndps_match <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                             treatment = "disaster_wo_operativo", refinement.method = "ps.weight", 
                             data = hrd, match.missing = T, 
                             covs.formula = ~  I(lag(subsetSECFORCESpc, 1:2)) + I(lag(lnpob, 1:2)) , 
                             qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                             exact.match.variables = c("gradorezagosoc2010num", "rural"),
                             lead = 0:0, forbid.treatment.reversal = FALSE)

fourndps_estim <- PanelEstimate(sets = fourndps_match, data = hrd)

summary(fourndps_estim)
fournd_ps<-as.data.frame(summary(fourndps_estim))


# 5. TWO LAG, ONE LEAD

fivendps_match <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                             treatment = "disaster_wo_operativo", refinement.method = "ps.weight", 
                             data = hrd, match.missing = T, 
                             covs.formula = ~  I(lag(subsetSECFORCESpc, 1:2)) + I(lag(lnpob, 1:2)) , 
                             qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                             exact.match.variables = c("gradorezagosoc2010num", "rural"),
                             lead = 0:1, forbid.treatment.reversal = FALSE)

fivendps_estim <- PanelEstimate(sets = fivendps_match, data = hrd)

summary(fivendps_estim)
fivend_ps<-as.data.frame(summary(fivendps_estim))



# 6. TWO LAG, TWO LEAD

sixndps_match <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                            treatment = "disaster_wo_operativo", refinement.method = "ps.weight", 
                            data = hrd, match.missing = T, 
                            covs.formula = ~ I(lag(subsetSECFORCESpc, 1:2)) + I(lag(lnpob, 1:2)) , 
                            qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                            exact.match.variables = c("gradorezagosoc2010num", "rural"),
                            lead = 0:2, forbid.treatment.reversal = FALSE)

sixndps_estim <- PanelEstimate(sets = sixndps_match, data = hrd)

summary(sixndps_estim)
sixnd_ps<-as.data.frame(summary(sixndps_estim))


rm(onendps_estim,onendps_match,twondps_estim,twondps_match,threendps_estim,threendps_match,
   fourndps_estim,fourndps_match,fivendps_estim,fivendps_match,sixndps_estim,sixndps_match)

##------------ CBPS WEIGHTS ####

# 1. ONE LAG, ZERO LEAD

onendcpbs_match <- PanelMatch(lag = 1, time.id = "year", unit.id = "inegi", 
                              treatment = "disaster_wo_operativo", refinement.method = "CBPS.weight", 
                              data = hrd, match.missing = T, 
                              covs.formula = ~ I(lag(subsetSECFORCESpc, 1:1)) + I(lag(lnpob, 1:1)) , 
                              size.match = 5, qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                              exact.match.variables = c("gradorezagosoc2010num", "rural"),
                              lead = 0:0, forbid.treatment.reversal = FALSE)

onendcbps_estim <- PanelEstimate(sets = onendcpbs_match, data = hrd)

summary(onendcbps_estim)
onend_cbps<-as.data.frame(summary(onendcbps_estim))

# 2. ONE LAG, ONE LEADS

twondcbps_match <- PanelMatch(lag = 1, time.id = "year", unit.id = "inegi", 
                              treatment = "disaster_wo_operativo", refinement.method = "CBPS.weight", 
                              data = hrd, match.missing = T, 
                              covs.formula = ~ I(lag(subsetSECFORCESpc, 1:1)) + I(lag(lnpob, 1:1)) , 
                              size.match = 5, qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                              exact.match.variables = c("gradorezagosoc2010num", "rural"),
                              lead = 0:1, forbid.treatment.reversal = FALSE)

twondcbps_estim <- PanelEstimate(sets = twondcbps_match, data = hrd)

summary(twondcbps_estim)
twond_cbps<-as.data.frame(summary(twondcbps_estim)) 

# 3. ONE LAG, TWO LEADS

threendcbps_match <- PanelMatch(lag = 1, time.id = "year", unit.id = "inegi", 
                                treatment = "disaster_wo_operativo", refinement.method = "CBPS.weight", 
                                data = hrd, match.missing = T, 
                                covs.formula = ~ I(lag(subsetSECFORCESpc, 1:1)) + I(lag(lnpob, 1:1)) , 
                                qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                                exact.match.variables = c("gradorezagosoc2010num", "rural"),
                                lead = 0:2, forbid.treatment.reversal = FALSE)

threendcbps_estim <- PanelEstimate(sets = threendcbps_match, data = hrd)

summary(threendcbps_estim)
threend_cbps<-as.data.frame(summary(threendcbps_estim))

# 4. TWO LAG, ZERO LEAD

fourndcbps_match <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                               treatment = "disaster_wo_operativo", refinement.method = "CBPS.weight", 
                               data = hrd, match.missing = T, 
                               covs.formula = ~ I(lag(subsetSECFORCESpc, 1:2)) + I(lag(lnpob, 1:2)) , 
                               qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                               exact.match.variables = c("gradorezagosoc2010num", "rural"),
                               lead = 0:0, forbid.treatment.reversal = FALSE)

fourndcbps_estim <- PanelEstimate(sets = fourndcbps_match, data = hrd)

summary(fourndcbps_estim)
fournd_cbps<-as.data.frame(summary(fourndcbps_estim))

# 5. TWO LAG, ONE LEAD

fivendcbps_match <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                               treatment = "disaster_wo_operativo", refinement.method = "CBPS.weight", 
                               data = hrd, match.missing = T, 
                               covs.formula = ~ I(lag(subsetSECFORCESpc, 1:2)) + I(lag(lnpob, 1:2)) , 
                               qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                               exact.match.variables = c("gradorezagosoc2010num", "rural"),
                               lead = 0:1, forbid.treatment.reversal = FALSE)

fivendcbps_estim <- PanelEstimate(sets = fivendcbps_match, data = hrd)

summary(fivendcbps_estim)
fivend_cbps<-as.data.frame(summary(fivendcbps_estim))


# 6. TWO LAG, TWO LEAD

sixndcbps_match <- PanelMatch(lag = 2, time.id = "year", unit.id = "inegi", 
                              treatment = "disaster_wo_operativo", refinement.method = "CBPS.weight", 
                              data = hrd, match.missing = T, 
                              covs.formula = ~  I(lag(subsetSECFORCESpc, 1:2)) + I(lag(lnpob, 1:2)) , 
                              qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                              exact.match.variables = c("gradorezagosoc2010num", "rural"),
                              lead = 0:2, forbid.treatment.reversal = FALSE)

sixndcbps_estim <- PanelEstimate(sets = sixndcbps_match, data = hrd)

summary(sixndcbps_estim)
sixnd_cbps<-as.data.frame(summary(sixndcbps_estim))


rm(onendcpbs_match,onendcbps_estim,twondcbps_estim,twondcbps_match,threendcbps_estim,threendcbps_match,
   fourndcbps_estim,fourndcbps_match,fivendcbps_estim,fivendcbps_match,sixndcbps_estim,sixndcbps_match)

##------------ FIGURE ####

# GENERATE NEW VARIABLES AND RENAME

onend_ps$Method<-c("PS Weights")
onend_cbps$Method<-c("CBPS Weights")
twond_ps$Method<-c("PS Weights")
twond_cbps$Method<-c("CBPS Weights")
threend_ps$Method<-c("PS Weights")
threend_cbps$Method<-c("CBPS Weights")
fournd_ps$Method<-c("PS Weights")
fournd_cbps$Method<-c("CBPS Weights")
fivend_ps$Method<-c("PS Weights")
fivend_cbps$Method<-c("CBPS Weights")
sixnd_ps$Method<-c("PS Weights")
sixnd_cbps$Method<-c("CBPS Weights")

onend_ps<-onend_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

onend_cbps<-onend_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

twond_ps<-twond_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

twond_cbps<-twond_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

threend_ps<-threend_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

threend_cbps<-threend_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

fournd_ps<-fournd_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

fournd_cbps<-fournd_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

fivend_ps<-fivend_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

fivend_cbps<-fivend_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

sixnd_ps<-sixnd_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )

sixnd_cbps<-sixnd_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5.,
    Lag=lag
  )


onend_ps$Time<-c("t+0")
onend_cbps$Time<-c("t+0")
twond_ps$Time<-c("t+0","t+1")
twond_cbps$Time<-c("t+0","t+1")
threend_ps$Time<-c("t+0","t+1","t+2")
threend_cbps$Time<-c("t+0","t+1","t+2")
fournd_ps$Time<-c("t+0")
fournd_cbps$Time<-c("t+0")
fivend_ps$Time<-c("t+0","t+1")
fivend_cbps$Time<-c("t+0","t+1")
sixnd_ps$Time<-c("t+0","t+1","t+2")
sixnd_cbps$Time<-c("t+0","t+1","t+2")

onend_ps$Model<-c("1 Lag 0 Lead")
onend_cbps$Model<-c("1 Lag 0 Lead")
twond_ps$Model<-c("1 Lag 1 Lead")
twond_cbps$Model<-c("1 Lag 1 Lead")
threend_ps$Model<-c("1 Lag 2 Lead")
threend_cbps$Model<-c("1 Lag 2 Lead")
fournd_ps$Model<-c("2 Lag 0 Lead")
fournd_cbps$Model<-c("2 Lag 0 Lead")
fivend_ps$Model<-c("2 Lag 1 Lead")
fivend_cbps$Model<-c("2 Lag 1 Lead")
sixnd_ps$Model<-c("2 Lag 2 Lead")
sixnd_cbps$Model<-c("2 Lag 2 Lead")

# APPEND DATA SETS


lagleadnd<-rbind(onend_cbps,onend_ps,twond_cbps,twond_ps,threend_cbps,threend_ps,
                 fournd_cbps,fournd_ps,fivend_cbps,fivend_ps,sixnd_cbps,sixnd_ps)

# CREATE FIGURE

theme_bw2 <- function(base_size = 13, base_family = "") {
  theme_grey(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = 8, colour = "black",  hjust = .5 , vjust=1),
      axis.text.y =       element_text(size = 11 , colour = "black", hjust = 0 , vjust=.5 ), # changes position of X axis text
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90),
      axis.title.x =      element_text(size =12),
      legend.position = "bottom",
      panel.background = element_rect(fill = "white"),
      strip.background = element_rect(colour="black", fill="white"),
      legend.title=element_blank(),
      legend.text=element_text(size=7),
      legend.key.size = unit(.7,"line")
    )
}

lagleadsnd_fig=ggplot(lagleadnd, aes(Time, Estimate, fill=Model)) + 
  geom_linerange(aes(ymin=CILow, ymax=CIHigh, colour=Model), size=1,  position = position_dodge(width = 0.5)) + 
  ylim(-2, 2) +
  theme_bw1() +
  geom_hline(yintercept=0, color="gray",  linetype=2) +
  facet_grid(.~Method)

plot(lagleadsnd_fig)  

ggsave("2_FIGURES/SA_FIG9A.png", plot=lagleadsnd_fig, width=8,height=4, units = "in", dpi = 300)

rm(onend_cbps,onend_ps,twond_cbps,twond_ps,threend_cbps,threend_ps,fournd_cbps,fournd_ps,
   fivend_cbps,fivend_ps,sixnd_cbps,sixnd_ps)

#### ------ STEP 12: SUPPLEMENTAL APPENDIX FIGURE 10A ####
##------------ PS WEIGHTS ####

# MATCHING
psmen_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                          treatment = "operativo", refinement.method = "ps.weight", 
                          data = hr, match.missing = T, 
                          covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(pobh15a44pc, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                          qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                          exact.match.variables = c("gradorezagosoc2010num", "rural"),
                          lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE
get_covariate_balance(psmen_match$att, hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "pobh15a44pc","turfwar_threeaf","samegovpres"), plot = F, ylim = c(-2,2))

# DIFF & DIFF
psmen_estim <- PanelEstimate(sets = psmen_match, data = hr)
summary(psmen_estim)
men_ps<-as.data.frame(summary(psmen_estim))

##------------ CBPS WEIGHTS ####

# MATCHING
cbpsmen_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                            treatment = "operativo", refinement.method = "CBPS.weight", 
                            data = hr, match.missing = T, 
                            covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(pobh15a44pc, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                            qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                            exact.match.variables = c("gradorezagosoc2010num", "rural"),
                            lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE
get_covariate_balance(cbpsmen_match$att, hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "pobh15a44pc","turfwar_threeaf","samegovpres"), plot = F, ylim = c(-2,2))

# DIFF & DIFF
cbpsmen_estim <- PanelEstimate(sets = cbpsmen_match, data = hr)
summary(cbpsmen_estim)
men_cbps<-as.data.frame(summary(cbpsmen_estim))

##------------ FIGURE ####

# GENERATE NEW VARIABLES AND RENAME

men_ps$Method<-c("PS weights")
men_cbps$Method<-c("CBPS weights")

men_ps<-men_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

men_cbps<-men_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

men_ps$Time<-c("t+0","t+1","t+2","t+3" )
men_cbps$Time<-c("t+0","t+1","t+2","t+3" )

# APPEND DATA SETS

men<-rbind(men_cbps,men_ps)

# CREATE FIGURE

men_fig=ggplot(men, aes(Time, Estimate)) + 
  geom_point() + 
  geom_linerange(aes(ymin=CILow, ymax=CIHigh), size=1, color="black") + 
  facet_grid(.~Method) + 
  ylim(-2, 2) +
  theme_bw1() +
  geom_hline(yintercept=0, color="gray",  linetype=2)

plot(men_fig)

ggsave("2_FIGURES/SA_FIG10A.png", plot=men_fig, width=6,height=4, units = "in", dpi = 300)

rm(men_cbps,men_ps,psmen_estim,psmen_match,cbpsmen_estim,cbpsmen_match)

### ------ STEP 13: SUPPLEMENTAL APPENDIX FIGURE 11A ####

hr_wosemar<-subset(hr, t3_semar_alone==0 & t6_sem_polfed==0)

##------------ PS WEIGHTS ####

# MATCHING
wosemarps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                              treatment = "operativo", refinement.method = "ps.weight", 
                              data = hr_wosemar, match.missing = T, 
                              covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                              qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                              exact.match.variables = c("gradorezagosoc2010num", "rural"),
                              lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE
get_covariate_balance(wosemarps_match$att, hr_wosemar, covariates = c("thom", "subsetSECFORCESpc", "lnpob","turfwar_threeaf","samegovpres"), plot = F, ylim = c(-2,2))


# DIFF & DIFF
wosemarps_estim <- PanelEstimate(sets = wosemarps_match, data = hr_wosemar)

summary(wosemarps_estim)

wosemar_ps<-as.data.frame(summary(wosemarps_estim))


##------------ CBPS WEIGHTS ####

# MATCHING
wosemarcbps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                                treatment = "operativo", refinement.method = "CBPS.weight", 
                                data = hr_wosemar, match.missing = T, 
                                covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                                qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                                exact.match.variables = c("gradorezagosoc2010num", "rural"),
                                lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE
get_covariate_balance(wosemarcbps_match$att, hr_wosemar, covariates = c("thom", "subsetSECFORCESpc", "lnpob","turfwar_threeaf","samegovpres"), plot = F, ylim = c(-2,2))


# DIFF & DIFF
wosemarcbps_estim <- PanelEstimate(sets = wosemarcbps_match, data = hr_wosemar)

summary(wosemarcbps_estim)

wosemar_cbps<-as.data.frame(summary(wosemarcbps_estim))

##------------ FIGURE ####

# GENERATE NEW VARIABLES AND RENAME

wosemar_ps$Method<-c("PS weights")
wosemar_cbps$Method<-c("CBPS weights")

wosemar_ps<-wosemar_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

wosemar_cbps<-wosemar_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

wosemar_ps$Time<-c("t+0","t+1","t+2","t+3" )
wosemar_cbps$Time<-c("t+0","t+1","t+2","t+3" )

# APPEND DATA SETS

wosemar<-rbind(wosemar_cbps,wosemar_ps)

# CREATE FIGURE

wosemar_fig=ggplot(wosemar, aes(Time, Estimate)) + 
  geom_point() + 
  geom_linerange(aes(ymin=CILow, ymax=CIHigh), size=1, color="black") + 
  facet_grid(.~Method) + 
  ylim(-2, 2) +
  theme_bw1() +
  geom_hline(yintercept=0, color="gray",  linetype=2)

plot(wosemar_fig)

ggsave("2_FIGURES/SA_FIG11A.png", plot=wosemar_fig, width=6,height=4, units = "in", dpi = 300)

rm(wosemar_cbps,wosemar_ps,wosemarcbps_estim,wosemarcbps_match,wosemarps_estim,wosemarps_match)



#### ------ STEP 14: SUPPLEMENTAL APPENDIX FIGURE 12A ####
##------------ PS WEIGHTS ####

# MATCHING
sedpfps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                            treatment = "operativo", refinement.method = "ps.weight", 
                            data = hr, match.missing = T, 
                            covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSEDPFpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                            qoi = "att" ,outcome.var = "subsetSEDPFpc",
                            exact.match.variables = c("gradorezagosoc2010num", "rural"),
                            lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE
get_covariate_balance(sedpfps_match$att, hr, covariates = c("thom", "subsetSEDPFpc", "lnpob","turfwar_threeaf","samegovpres"), plot = F, ylim = c(-2,2))


# DIFF & DIFF
sedpfps_estim <- PanelEstimate(sets = sedpfps_match, 
                               data = hr)
summary(sedpfps_estim)

sedpf_ps<-as.data.frame(summary(sedpfps_estim))


##------------ CBPS WEIGHTS ####

# MATCHING
sedpfcbps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                              treatment = "operativo", refinement.method = "CBPS.weight", 
                              data = hr, match.missing = T, 
                              covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSEDPFpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                              qoi = "att" ,outcome.var = "subsetSEDPFpc",
                              exact.match.variables = c("gradorezagosoc2010num", "rural"),
                              lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE
get_covariate_balance(sedpfcbps_match$att, hr, covariates = c("thom", "subsetSEDPFpc", "lnpob","turfwar_threeaf","samegovpres"), plot = F, ylim = c(-2,2))


# DIFF & DIFF
sedpfcbps_estim <- PanelEstimate(sets = sedpfcbps_match, data = hr)

summary(sedpfcbps_estim)

sedpf_cbps<-as.data.frame(summary(sedpfcbps_estim))


##------------ FIGURE ####

# GENERATE NEW VARIABLES AND RENAME

sedpf_ps$Method<-c("PS weights")
sedpf_cbps$Method<-c("CBPS weights")

sedpf_ps<-sedpf_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

sedpf_cbps<-sedpf_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

sedpf_ps$Time<-c("t+0","t+1","t+2","t+3" )
sedpf_cbps$Time<-c("t+0","t+1","t+2","t+3" )

# APPEND DATA SETS

sedpf<-rbind(sedpf_cbps,sedpf_ps)

# CREATE FIGURE

sedpf_fig=ggplot(sedpf, aes(Time, Estimate)) + 
  geom_point() + 
  geom_linerange(aes(ymin=CILow, ymax=CIHigh), size=1, color="black") + 
  facet_grid(.~Method) + 
  ylim(-2, 2) +
  theme_bw1() +
  geom_hline(yintercept=0, color="gray",  linetype=2)

plot(sedpf_fig)

ggsave("2_FIGURES/SA_FIG12A.png", plot=sedpf_fig, width=6,height=4, units = "in", dpi = 300)

rm(sedpf_cbps,sedpf_ps,sedpfcbps_estim,sedpfcbps_match,sedpfps_estim,sedpfps_match)

#### ------ STEP 15: SUPPLEMENTAL APPENDIX FIGURE 13A #### 
##------------ PS WEIGHTS ####

# MATCHING  

copmayps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                             treatment = "operativo", refinement.method = "ps.weight", 
                             data = hr, match.missing = T, 
                             covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samepartypres, 1:3)), 
                             qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                             exact.match.variables = c("gradorezagosoc2010num", "rural"),
                             lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE

get_covariate_balance(copmayps_match$att, hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "turfwar_threeaf", "samepartypres"), plot = F, ylim = c(-2,2))

# DIFF & DIFF

copmayorps_estim <- PanelEstimate(sets = copmayps_match, data = hr)

summary(copmayorps_estim)

copmayor_ps<-as.data.frame(summary(copmayorps_estim))

##------------ CBPS WEIGHTS ####


# MATCHING 

copmayorcbps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                                 treatment = "operativo", refinement.method = "CBPS.weight", 
                                 data = hr, match.missing = T, 
                                 covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samepartypres, 1:3)), 
                                 qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                                 exact.match.variables = c("gradorezagosoc2010num", "rural"),
                                 lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE

get_covariate_balance(copmayorcbps_match$att, hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "turfwar_threeaf", "samepartypres"), plot = F, ylim = c(-2,2))

# DIFF & DIFF  

copmayorcbps_estim <- PanelEstimate(sets = copmayorcbps_match, data = hr)

summary(copmayorcbps_estim)

copmayor_cbps<-as.data.frame(summary(copmayorcbps_estim))

##------------ FIGURE ####

# GENERATE NEW VARIABLES AND RENAME

copmayor_ps$Method<-c("PS Weights")
copmayor_cbps$Method<-c("CBPS Weights")

copmayor_ps<-copmayor_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

copmayor_cbps<-copmayor_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

copmayor_ps$Time<-c("t+0","t+1","t+2","t+3" )
copmayor_cbps$Time<-c("t+0","t+1","t+2","t+3" )

# APPEND DATA SETS

copmayor<-rbind(copmayor_cbps,copmayor_ps)

# CREATE FIGURE

copmayor_fig=ggplot(copmayor, aes(Time, Estimate)) + 
  geom_point() + 
  geom_linerange(aes(ymin=CILow, ymax=CIHigh), size=1, color="black") + 
  facet_grid(.~Method) + 
  ylim(-2, 2) +
  theme_bw1() +
  geom_hline(yintercept=0, color="gray",  linetype=2)

plot(copmayor_fig)

ggsave("2_FIGURES/SA_FIG13A.png", plot=copmayor_fig, width=6,height=4, units = "in", dpi = 300)

rm(copmayor_cbps,copmayor_ps,copmayorcbps_estim,copmayorcbps_match,copmayorps_estim,copmayps_match)

#### ------ STEP 16: SUPPLEMENTAL APPENDIX FIGURE 14A ####
##------------ PS WEIGHTS ####

# MATCHING  

subsubps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                             treatment = "operativo", refinement.method = "ps.weight", 
                             data = hr, match.missing = T, 
                             covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsubsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                             qoi = "att" ,outcome.var = "subsubsetSECFORCESpc",
                             exact.match.variables = c("gradorezagosoc2010num", "rural"),
                             lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE

get_covariate_balance(subsubps_match$att, hr, covariates = c("thom", "subsubsetSECFORCESpc", "lnpob","turfwar_threeaf","samegovpres"), plot = F, ylim = c(-2,2))

# DIFF & DIFF

subsubps_estim <- PanelEstimate(sets = subsubps_match, data = hr)

summary(subsubps_estim)

subsub_ps<-as.data.frame(summary(subsubps_estim))

##------------ CBPS WEIGHTS ####


# MATCHING 

subsubcbps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                               treatment = "operativo", refinement.method = "CBPS.weight", 
                               data = hr, match.missing = T, 
                               covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsubsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                               qoi = "att" ,outcome.var = "subsubsetSECFORCESpc",
                               exact.match.variables = c("gradorezagosoc2010num", "rural"),
                               lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE

get_covariate_balance(subsubcbps_match$att, hr, covariates = c("thom", "subsubsetSECFORCESpc", "lnpob","turfwar_threeaf","samegovpres"), plot = F, ylim = c(-2,2))

# DIFF & DIFF  

subsubcbps_estim <- PanelEstimate(sets = subsubcbps_match, data = hr)

summary(subsubcbps_estim)

subsub_cbps<-as.data.frame(summary(subsubcbps_estim))

##------------ FIGURE ####

# GENERATE NEW VARIABLES AND RENAME

subsub_ps$Method<-c("PS Weights")
subsub_cbps$Method<-c("CBPS Weights")

subsub_ps<-subsub_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

subsub_cbps<-subsub_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

subsub_ps$Time<-c("t+0","t+1","t+2","t+3" )
subsub_cbps$Time<-c("t+0","t+1","t+2","t+3" )

# APPEND DATA SETS

subsub<-rbind(subsub_cbps,subsub_ps)

# CREATE FIGURE

subsub_fig=ggplot(subsub, aes(Time, Estimate)) + 
  geom_point() + 
  geom_linerange(aes(ymin=CILow, ymax=CIHigh), size=1, color="black") + 
  facet_grid(.~Method) + 
  ylim(-2, 2) +
  theme_bw1() +
  geom_hline(yintercept=0, color="gray",  linetype=2)

plot(subsub_fig)

ggsave("2_FIGURES/SA_FIG14A.png", plot=subsub_fig, width=6,height=4, units = "in", dpi = 300)

rm(subsubcbps_estim,subsubcbps_match,subsubps_estim,subsubps_match)

#### ------ STEP 17: SUPPLEMENTAL APPENDIX FIGURE 15A #### 
##------------ PS WEIGHTS ####

# MATCHING  

twosdps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                            treatment = "operativo", refinement.method = "ps.weight", 
                            data = hr, match.missing = T, 
                            covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_af, 1:3)) + I(lag(samegovpres, 1:3)), 
                            qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                            exact.match.variables = c("gradorezagosoc2010num", "rural"),
                            lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE

get_covariate_balance(twosdps_match$att, hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "turfwar_af", "samegovpres"), plot = F, ylim = c(-2,2))

# DIFF & DIFF

twosdps_estim <- PanelEstimate(sets = twosdps_match, data = hr)

summary(twosdps_estim)

twosd_ps<-as.data.frame(summary(twosdps_estim))

##------------ CBPS WEIGHTS ####


# MATCHING 

twosdcbps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                              treatment = "operativo", refinement.method = "CBPS.weight", 
                              data = hr, match.missing = T, 
                              covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSECFORCESpc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_af, 1:3)) + I(lag(samegovpres, 1:3)), 
                              qoi = "att" ,outcome.var = "subsetSECFORCESpc",
                              exact.match.variables = c("gradorezagosoc2010num", "rural"),
                              lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE

get_covariate_balance(twosdcbps_match$att, hr, covariates = c("thom", "subsetSECFORCESpc", "lnpob", "turfwar_af", "samegovpres"), plot = F, ylim = c(-2,2))

# DIFF & DIFF  

twosdcbps_estim <- PanelEstimate(sets = twosdcbps_match, data = hr)

summary(twosdcbps_estim)

twosd_cbps<-as.data.frame(summary(twosdcbps_estim))

##------------ FIGURE ####

# GENERATE NEW VARIABLES AND RENAME

twosd_ps$Method<-c("PS weights")
twosd_cbps$Method<-c("CBPS weights")

twosd_ps<-twosd_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

twosd_cbps<-twosd_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

twosd_ps$Time<-c("t+0","t+1","t+2","t+3" )
twosd_cbps$Time<-c("t+0","t+1","t+2","t+3" )

# APPEND DATA SETS

twosd<-rbind(twosd_cbps,twosd_ps)

# CREATE FIGURE

twosd_fig=ggplot(twosd, aes(Time, Estimate)) + 
  geom_point() + 
  geom_linerange(aes(ymin=CILow, ymax=CIHigh), size=1, color="black") + 
  facet_grid(.~Method) + 
  ylim(-2, 2) +
  theme_bw1() +
  geom_hline(yintercept=0, color="gray",  linetype=2)

plot(twosd_fig)

ggsave("2_FIGURES/SA_FIG15A.png", plot=twosd_fig, width=6,height=4, units = "in", dpi = 300)

rm(twosd_cbps,twosd_ps,twosdcbps_estim,twosdcbps_match,twosdps_estim,twosdps_match)

#### ------ STEP 18: SUPPLEMENTAL APPENDIX FIGURE 16A #### 
##------------ PS WEIGHTS ####

# MATCHING  

mainsaps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                           treatment = "operativo", refinement.method = "ps.weight", 
                           data = hr, match.missing = T, 
                           covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSEDENApc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                           qoi = "att" ,outcome.var = "subsetSEDENApc",
                           exact.match.variables = c("gradorezagosoc2010num", "rural"),
                           lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE

get_covariate_balance(mainsaps_match$att, hr, covariates = c("thom", "subsetSEDENApc", "lnpob", "turfwar_threeaf", "samegovpres"), plot = F, ylim = c(-2,2))

# DIFF & DIFF

mainsaps_estim <- PanelEstimate(sets = mainsaps_match, data = hr)

summary(mainsaps_estim)

mainsa_ps<-as.data.frame(summary(mainsaps_estim))

##------------ CBPS WEIGHTS ####


# MATCHING 

mainsacbps_match <- PanelMatch(lag = 3, time.id = "year", unit.id = "inegi", 
                             treatment = "operativo", refinement.method = "CBPS.weight", 
                             data = hr, match.missing = T, 
                             covs.formula = ~ I(lag(thom, 1:3))  + I(lag(subsetSEDENApc, 1:3)) + I(lag(lnpob, 1:3)) + I(lag(turfwar_threeaf, 1:3)) + I(lag(samegovpres, 1:3)), 
                             qoi = "att" ,outcome.var = "subsetSEDENApc",
                             exact.match.variables = c("gradorezagosoc2010num", "rural"),
                             lead = 0:3, forbid.treatment.reversal = FALSE)

# GET COVARIATE BALANCE

get_covariate_balance(mainsacbps_match$att, hr, covariates = c("thom", "subsetSEDENApc", "lnpob", "turfwar_threeaf", "samegovpres"), plot = F, ylim = c(-2,2))

# DIFF & DIFF  

mainsacbps_estim <- PanelEstimate(sets = mainsacbps_match, data = hr)

summary(mainsacbps_estim)

mainsa_cbps<-as.data.frame(summary(mainsacbps_estim))

##------------ FIGURE ####

# GENERATE NEW VARIABLES AND RENAME

mainsa_ps$Method<-c("PS weights")
mainsa_cbps$Method<-c("CBPS weights")

mainsa_ps<-mainsa_ps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

mainsa_cbps<-mainsa_cbps %>% 
  rename(
    Estimate=summary.estimate,
    StError=summary.std.error,
    CILow=summary.2.5.,
    CIHigh=summary.97.5. 
  )

mainsa_ps$Time<-c("t+0","t+1","t+2","t+3" )
mainsa_cbps$Time<-c("t+0","t+1","t+2","t+3" )

# APPEND DATA SETS

mainsa<-rbind(mainsa_cbps,mainsa_ps)

# CREATE FIGURE

mainsa_fig=ggplot(mainsa, aes(Time, Estimate)) + 
  geom_point() + 
  geom_linerange(aes(ymin=CILow, ymax=CIHigh), size=1, color="black") + 
  facet_grid(.~Method) + 
  ylim(-2, 2) +
  theme_bw1() +
  geom_hline(yintercept=0, color="gray",  linetype=2)

plot(mainsa_fig)

ggsave("2_FIGURES/SA_FIG16A.png", plot=mainsa_fig, width=6,height=4, units = "in", dpi = 300)

rm(mainsa_cbps,mainsa_ps,mainsacbps_estim,mainsacbps_match,mainsaps_estim,mainsaps_match)

